

! cut-off functions for robust estimation of statistical parameters

module cutoff

  integer, parameter, private :: dp = selected_real_kind(15)
  real(dp),parameter,private :: huber_a = 1.345
  real(dp),parameter,private :: hampel_a = 1.7, hampel_b = 3.4, hampel_c = 8.5
  real(dp),parameter,private :: Andrews_a = 2.1, Pi = 3.1415929
  real(dp),parameter,private :: Tukey_a = 6.0

contains 

! -------------------------------------------------------------
! 
! Huber's function

function huber(x)

  real(dp) :: x,huber
 
  if( abs(x) < huber_a )then
     huber = x
  elseif( abs(x) >= huber_a )then
     huber = sign(huber_a,x)
  else
     huber = sign(huber_a,x)
!     write(*,*) x,huber_a   ! x = NaN
!     stop '1'
  endif

end function huber

function dhuber(x)

  real(dp) :: x,dhuber
  
  if( abs(x) < huber_a )then ! absolute value of derivation of x
      dhuber = 1.0_dp
  elseif( abs(x) >= huber_a )then
     dhuber = 0.0_dp
  else
     dhuber = 0.0_dp
!     write(*,*) x,huber_a
!     stop '2'
  endif

end function dhuber
 
!-----------------------------------------------------------------------
!
! Hampel's function

function hampel(x)

  real(dp) :: x,hampel
 
  if( abs(x) < hampel_a )then
     hampel = x
  elseif( hampel_a <= abs(x) .and. abs(x) < hampel_b )then
     hampel = sign(hampel_a,x)
  elseif( hampel_b <= abs(x) .and. abs(x) < hampel_c )then
     hampel = hampel_a*(x - sign(hampel_c,x))/(hampel_b - hampel_c)
  elseif(  abs(x) >= hampel_c )then
     hampel = 0.0_dp
  else
     hampel = 0.0_dp
!     write(*,*) x,hampel_a
!     stop '3'
  endif

end function hampel

function dhampel(x)

  real(dp) :: x,dhampel
  
  if( abs(x) < hampel_a )then    ! absolute value of derivation of x
      dhampel = 1.0_dp
  elseif( hampel_b <= abs(x) .and. abs(x) < hampel_c )then
     dhampel = - sign(hampel_a/(hampel_b - hampel_c),x)
  else
     dhampel = 0.0_dp
  endif

end function dhampel

!-----------------------------------------------------------------------
!
! Andrews's function

function Andrews(x)

  real(dp) :: x,Andrews
 
  if( abs(x) < Andrews_a*Pi )then
     Andrews = sin(x/Andrews_a)
  else
     Andrews = 0.0_dp
  endif

end function Andrews

function dAndrews(x)

  real(dp) :: x,dAndrews
  
  if( abs(x) < Andrews_a*Pi )then ! absolute value of derivation of x
      dAndrews = cos(x/Andrews_a)/Andrews_a
  else
     dAndrews = 0.0_dp
  endif

end function dAndrews

!-----------------------------------------------------------------------
!
! Tukey's function

function Tukey(x)

  real(dp) :: x,Tukey
 
  if( abs(x) < Tukey_a )then
     Tukey = x*(1.0_dp - (x/Tukey_a)**2)**2
  else
     Tukey = 0.0_dp
  endif

end function Tukey

function dTukey(x)

  real(dp) :: x,dTukey
  
  if( abs(x) < Tukey_a )then ! absolute value of derivation of x
      dTukey = (1.0_dp - (x/Tukey_a)**2)*(1.0_dp - 5.0_dp*(x/Tukey_a)**2)
  else
     dTukey = 0.0_dp
  endif

end function dTukey

end module cutoff


